A chromatic feature detector in the retina signals visual context changes

The retina transforms patterns of light into visual feature representations supporting behaviour. These representations are distributed across various types of retinal ganglion cells (RGCs), whose spatial and temporal tuning properties have been studied extensively in many model organisms, including the mouse. However, it has been difficult to link the potentially nonlinear retinal transformations of natural visual inputs to specific ethological purposes. Here, we discover a nonlinear selectivity to chromatic contrast in an RGC type that allows the detection of changes in visual context. We trained a convolutional neural network (CNN) model on large-scale functional recordings of RGC responses to natural mouse movies, and then used this model to search in silico for stimuli that maximally excite distinct types of RGCs. This procedure predicted centre colour opponency in transient suppressed-by-contrast (tSbC) RGCs, a cell type whose function is being debated. We confirmed experimentally that these cells indeed responded very selectively to Green-OFF, UV-ON contrasts. This type of chromatic contrast was characteristic of transitions from ground to sky in the visual scene, as might be elicited by head or eye movements across the horizon. Because tSbC cells performed best among all RGC types at reliably detecting these transitions, we suggest a role for this RGC type in providing contextual information (i.e. sky or ground) necessary for the selection of appropriate behavioural responses to other stimuli, such as looming objects. Our work showcases how a combination of experiments with natural stimuli and computational modelling allows discovering novel types of stimulus selectivity and identifying their potential ethological relevance.


Introduction
Sensory systems evolved to generate representations of an animal's natural environment useful for survival and procreation (1).These environments are complex and highdimensional, and different features are relevant for different species (reviewed in (2)).As a consequence, the representations are adapted to an animal's needs: features of the world relevant for the animal are represented with enhanced preci-sion, whereas less important features are discarded.Sensory (24)).In recent years, convolutional neural network (CNN) models have become the state-of-the-art approach for predictive modelling of visual processing, both in the retina (25-28), as well as in higher visual areas (29-31).In the cortex, two recent studies took the CNN modelling approach further, beyond response prediction, by probing the networks for stimuli that would maximally excite the modelled neurons (32,33).The resulting maximally exciting inputs (MEIs) were more complex and diverse than expected based on previous results obtained with synthetic stimuli and linear methods.Leveraging the power of this approach, another study highlighted the ethological relevance of colour by uncovering a state-dependent shift in chromatic preference of mouse V1 neurons, a shift that could facilitate the detection of aerial predators against a UV-bright sky (34).
Here, we combined the power of CNN-based modelling with large-scale recordings from RGCs to investigate colour processing in the mouse retina under natural stimulus conditions.Since mouse photoreceptors are sensitive to green and UV light (35), we recorded RGC responses to stimuli capturing the chromatic composition of natural mouse environments in these two chromatic channels.A model-guided search for MEIs in chromatic stimulus space predicted a novel type of chromatic tuning in transient Suppressed-by-Contrast (tSbC) RGCs, a type whose function is being debated (36)(37)(38).
A detailed in-silico characterisation followed up by experimental validation ex-vivo confirmed this cell type's pronounced and unique selectivity for dynamic full-field changes from green-dominated to UV-dominated scenes, a type of visual input that matches the scene statistics of transitions across the horizon (5,39,40).We therefore suggest a role for tSbC RGCs in detecting behaviourally relevant changes in visual context, such as a transitions from ground  1), performs convolutions with 3D space-time separable filters (2) followed by a nonlinear activation function (3) in two consecutive layers (2)(3)(4) within its core, and feeds the output of its core into a per-neuron readout.For each RGC, the readout convolves the feature maps with a learned RF modelled as a 2D Gaussian (5), and finally feeds a weighted sum of the resulting vector through a softplus nonlinearity (6) to yield the firing rate prediction for that RGC (7)  The responses elicited by the synthetic stimuli and the natural movie were diverse, displaying ON (Figure 1d, rows 4-9), ON-OFF (row 3) and OFF (rows 1 and 2), as well as sustained and transient characteristics (e.g., rows 8 and 4, respectively).Some responses were suppressed by temporal contrast (generally, rows 10, 11; at high contrast and frequency, row 9).A total of 6,984 GCL cells passed our response quality criteria (see Methods); 3,527 cells could be assigned to one of 32 previously characterised functional RGC groups (6) based on their responses to the chirp and moving bar stimuli using our recently developed classifier (Figure 1e; Figure 1-figure supplement 1Ia) (5).Cells assigned to any of groups 33-46 were considered displaced amacrine cells and were not analysed in this study (for detailed filtering pipeline, see Figure 1-figure supplement 1Ic).

CNN model captures diverse tuning of RGC groups and predicts MEIs.
We trained a CNN model on the RGCs' movie responses (Figure 2a) and evaluated model performance as the correlation between predicted and trial-     We presented these MEIs on a regularly spaced 5 × 5 333 grid to achieve approximate centring of stimuli on RGC RFs 334 in the recording field (Figure 4b,c).For these recordings,

335
we fit models whose readout parameters allowed us to es- these experiments with precise positioning of stimuli on the cells' RFs elicited the same responses as the 2P experiments confirms the validity of the grid-approach for stimulus presentation used in the latter.
Chromatic contrast selectivity derives from a nonlinear transformation of stimulus space.Next, we asked whether G 28 (tSbC) RGC's selectivity is a linear feature, as could be achieved by two linear filters with opposite signs for the two colour channels, or whether it is a nonlinear feature.To address this question, we tested whether an LN model (implemented using convolutions; see Methods) could recover the chromatic selectivity of G 28 by predicting MEIs using the LN model (Figure 6).We found that the LN model predicted colour-opponent MEIs for only 9 out of 36 (25%) G 28 RGCs (nonlinear CNN: 24 out of 36 (66%) colour-opponent MEIs; Figure 6a-c).This finding argues against the possibility that G 28 's colour opponency can be explained on the computational level by two opposite-sign linear filters operating on the two colour channels, which could be recovered by a LN model.Instead, it suggests the presence of a nonlinear dependency between chromatic contrast (of the stimulus) and chromatic selectivity (of the cell).
In other words, G 28 RGCs process stimuli differently depending on their chromatic contrast, a nonlinear feature that cannot be accurately captured by a LN model that makes a single estimate of the linear filter for the whole stimulus space.
To understand the nature of this dependency, we expanded the estimate of the model RGCs' tuning to colour contrast around the maximum (the MEI).We did this by mapping the model neurons' response and its gradient in 2D 435 chromatic contrast space (Figure 6c).This analysis revealed  To test if G 28 (tSbC) RGCs respond to such a stimulus, we used the transitions between movie clips (inter-clip transitions; cf. Figure 1b) as a proxy for the type of visual input elicited by head or eye movements: ground-to-ground and sky-to-sky transitions for horizontal movements without change in visual context, and ground-to-sky and sky-toground transitions for vertical movements with a change in visual context.We then calculated the contrast of these transitions in the green and UV channel and mapped them to the chromatic contrast stimulus space (Figure 7a).We found that ground-to-ground and sky-to-sky transitions were distributed along the diagonal, whereas the two transitions re-   Confidence is defined as the probability assigned to the predicted class by the random forest classifier (see ( 81)), and the threshold was set at ≥ 0.25.
The fourth step made sure only cells with a sufficient model prediction performance, defined as an average singletrial test set correlation of C(r (n) , r All cells passing steps 1-3 were included in the horizon detection analysis (Figure 7); all cells passing steps 1-4 were included in the MEI analysis (Figure 3); the "red" cells passing steps 1-4 were included in the MEI validation analysis (Figure 4).In the process of analysing MEIs, we fitted DoGs For electrophysiology experiments, stimuli were presented using a digital projector (DPM-FE4500MKII, EKB Thus, (2) is the temporal kernel parameterisation, that allows the 1046 model to learn a shared temporal filter that is made faster 1047 or slower for each specific scan (86).where N is the number of neurons, r (n) is the mea- white noise, and then iteratively updated x (n) i according to the gradient of the model neuron's response: where λ = 10 was the learning rate.The optimisation was performed using Stochastic Gradient Descent (SGD), and was subject to a norm and a range constraint.The norm constraint was applied jointly across both channels and ensured that the L2 norm of each MEI did not exceed a fixed budget b of 30.The norm-constrained MEI x(n) i was calculated at each iteration as The range constraint was defined and applied for each Analysing MEIs.We analysed MEIs to quantify their spatial, temporal, and chromatic properties.1159 and and likewise for G s .We initialised (µ x , µ y ) in the fol-  Butterworth filter with cutoff frequency 10 Hz.We then estimated the mean frequency of the temporal component by calculating an average of the frequency components, each weighted with its relative power.
Contrast The contrast of the MEIs in the two channels, ) for c ∈ [green, UV], was defined as the difference between the mean value within the centre mask m at the two last peaks of the temporal component of the MEI in the UV channel at time points t 2 and t 1 : where denotes the element-wise multiplication of the MEI and the binary mask.(see Figure 3f).The peaks were found with the function scipy.signal.find_peaks,and the peaks found for the UV channel were used to calculate contrast both in the green and the UV channel.

Validating MEIs experimentally.
Generating MEI stimuli.To test experimentally whether the model correctly predicts which stimuli would maximally excite RGCs of different RGC groups, we performed a new set of experiments (numbers indicated in red in Figure 1-figure supplement 1Ic), where we complemented our stimulus set with MEI stimuli.For the MEI stimuli, we selected 11 RGCs, chosen to span the responses space and to represent both well-described and poorly understood RGC groups, for which we generated MEIs at different positions on a 5 × 5 grid (spanning 110µm in vertical and horizontal direction).

Figure 1 .
Figure 1.Mouse RGCs display diverse responses to a natural movie stimulus (a) Illustration of a flat-mounted retina, with recording fields (white circles) and stimulus area centred on the red recording field indicated (cross marks optic disc; d, dorsal; v, ventral; t, temporal; n, nasal).(b) Natural movie stimulus structure (top) and example frames (bottom).The stimulus consisted of 5-s clips taken from UV-green footage recorded outside (5), with 3 repeats of a 5-clip test sequence (highlighted in grey) and a 108-clip training sequence (see Methods).(c) Representative recording field (bottom; marked by red square in (a)) showing somata of ganglion cell layer (GCL) cells loaded with Ca 2+ indicator OGB-1.(d) Ca 2+ responses of exemplary RGCs (indicated by circles in (c)) to chirp (left), moving bar (centre), and natural movie (right) stimulus.(e) Same recording field as in (c) but with cells colour-coded by functional RGC group (left; see Methods and (6)) and group responses (coloured, mean ± SD across cells; trace of example cells in (d) overlaid in black).

Figure 2 .
Figure 2. CNN model captures diverse tuning of RGC groups and predicts MEIs (a) Illustration of the CNN model and its output.The model takes natural movie clips as input (1), performs convolutions with 3D space-time separable filters (2) followed by a nonlinear activation function (3) in two consecutive layers (2-4) within its core, and feeds the output of its core into a per-neuron readout.For each RGC, the readout convolves the feature maps with a learned RF modelled as a 2D Gaussian (5), and finally feeds a weighted sum of the resulting vector through a softplus nonlinearity (6) to yield the firing rate prediction for that RGC (7).Numbers indicate averaged single-trial test set correlation between predicted (red) and recorded (black) responses.(b) Test set correlation between model prediction and neural response (averaged across three repetitions) as a function of response reliability (see Methods) for N=3,527 RGCs.Coloured dots correspond to example cells shown in Figure 1c-e.Dots in darker grey correspond to the N=1,947 RGCs that passed the model test correlation and movie response quality criterion (see Methods and Figure 1-figure supplement 1Ic).(c) Test set correlation (as in (b)) of model vs. test set correlation of a linearised version of the CNN model (for details, see Methods).Coloured dots correspond to RGC groups 1-32 (6).Dark and light grey dots as in (b).(d) Illustration of model-guided search for maximally exciting inputs (MEIs).The trained model captures neural tuning to stimulus features (far left; heat map illustrates "landscape´´of neural tuning to stimulus features).Starting from a randomly initialised input (2nd from left; a 3D tensor in space and time; only one colour channel illustrated here), the model follows the gradient along the tuning surface (far left) to iteratively update the input until it arrives at the stimulus (bottom right) that maximises the model neuron's activation within an optimisation time window (0.66 s, grey box, top right).
Figure 2. CNN model captures diverse tuning of RGC groups and predicts MEIs (a) Illustration of the CNN model and its output.The model takes natural movie clips as input (1), performs convolutions with 3D space-time separable filters (2) followed by a nonlinear activation function (3) in two consecutive layers (2-4) within its core, and feeds the output of its core into a per-neuron readout.For each RGC, the readout convolves the feature maps with a learned RF modelled as a 2D Gaussian (5), and finally feeds a weighted sum of the resulting vector through a softplus nonlinearity (6) to yield the firing rate prediction for that RGC (7).Numbers indicate averaged single-trial test set correlation between predicted (red) and recorded (black) responses.(b) Test set correlation between model prediction and neural response (averaged across three repetitions) as a function of response reliability (see Methods) for N=3,527 RGCs.Coloured dots correspond to example cells shown in Figure 1c-e.Dots in darker grey correspond to the N=1,947 RGCs that passed the model test correlation and movie response quality criterion (see Methods and Figure 1-figure supplement 1Ic).(c) Test set correlation (as in (b)) of model vs. test set correlation of a linearised version of the CNN model (for details, see Methods).Coloured dots correspond to RGC groups 1-32 (6).Dark and light grey dots as in (b).(d) Illustration of model-guided search for maximally exciting inputs (MEIs).The trained model captures neural tuning to stimulus features (far left; heat map illustrates "landscape´´of neural tuning to stimulus features).Starting from a randomly initialised input (2nd from left; a 3D tensor in space and time; only one colour channel illustrated here), the model follows the gradient along the tuning surface (far left) to iteratively update the input until it arrives at the stimulus (bottom right) that maximises the model neuron's activation within an optimisation time window (0.66 s, grey box, top right).
presentation, thereby allowing to assess the reliability of neuronal responses across the recording (Figure1b, top).

Figure 3 .
Figure 3. Spatial, temporal and chromatic properties of MEIs differ between RGC groups (a) Spatial component of three example MEIs for green (top), UV (middle) and overlay (bottom).Solid and dashed circles indicate MEI centre and surround fit, respectively.For display, spatial components s in the two channels were re-scaled to a similar range and displayed on a common grey-scale map ranging from black for −max(|s|) to white for max(|s|), i.e. symmetric about 0 (grey).(b) Spatio-temporal (y-t) plot for the three example MEIs (from (a)) at a central vertical slice for green (top), UV (middle) and overlay (bottom).Grey-scale map analogous to (a).(c) Trajectories through colour space over time for the centre of the three MEIs.Trajectories start at the origin (grey level); direction of progress indicated by arrow heads.Bottom right: Bounding boxes of the respective trajectory plots.(d) Calculation of MEI centre size, defined as σx + σy, with σx and σy the s.d. in horizontal and vertical direction, respectively, of the DoG fit to the MEI.(e) Calculation of MEI temporal frequency: Temporal components are transformed using Fast Fourier Transform, and MEI frequency is defined as the amplitude-weighted average frequency of the Fourier-transformed temporal component.(f) Calculation of centre contrast, which is defined as the difference in intensity at the last two peaks (indicated by t1 and t2, respectively, in (c)).For the example cell (orange markers and lines), green intensity decreases, resulting in OFF contrast, and UV intensity increases, resulting in ON contrast.(g) Distribution of green and UV MEI centre sizes across N=1,613 cells (example MEIs from (a-c) indicated by arrows; symbolsas shown on top of (a)).95% of MEIs were within an angle of ± 8°of the diagonal (solid and dashed lines); MEIs outside of this range are coloured by cell type.(h) As (g) but for distribution of green and UV MEI temporal frequency.95% of MEIs were within an angle of ± 11.4°of the diagonal (solid and dashed lines).(i) As (g) but for distribution of green and UV MEI centre contrast.MEI contrast is shifted away from the diagonal (dashed line) towards UV by an angle of 33.2°due to the dominance of UV-sensitive S-opsin in the ventral retina.MEIs at an angle > 45°occupy the upper left, colour-opponent (UV ON -green OFF ) quadrant.(j, k) Fraction of MEIs per cell type that lie outside the angle about the diagonal containing 95% of MEIs for centre size and temporal frequency.Broad RGC response types indicated as in (6).(l) Fraction of MEIs per cell type in the upper-left, colour-opponent quadrant for contrast.

Figure 4 .
Figure 4. Experiments confirm MEIs predicted by model (a) MEIs shown during the experiment, with green and UV spatial components (top two rows), as well as green and UV temporal components (third row) and a spatio-temporal visualisation (fourth row).For display, spatial components s in the two channels were re-scaled to a similar range and displayed on a common grey-scale map ranging from black for −max(|s|) to white for max(|s|), i.e. symmetric about 0 (grey).Relative amplitudes of UV and green are shown in the temporal components.(b) Illustration of spatial layout of MEI experiment.White circles represent 5 × 5 grid of positions where MEIs were shown; red shading shows an example RF estimate of a recorded G32 RGC, with black dot indicating the RF centre position (Methods).(c) Responses of example RGC from (b) to the 11 different MEI stimuli at 25 different positions.(d) Recorded (top, r (n) ) and predicted (bottom, r(n) ) responses to the 11 different MEIs for example RGC n from (b, c).Left: responses are averaged across the indicated dimensions x, y (different MEI locations); black bar indicates MEI stimulus duration (from 0 to 1.66 s), grey rectangle marks optimisation time window (from 1 to 1.66 s).Right: Response to different MEIs, additionally averaged across time (t; within optimisation time window).(e,f) Same as in (d), but additionally averaged across all RGCs (n) of G5 (N=6) (e) and of G28 (N=12) (f).Error bars show SD across cells.(g) Confusion matrix, each row showing the z-scored response magnitude of one RGC group (averaged across all RGCs of that group) to the MEIs in (a).Confusion matrix for recorded cells (top; "Data") and for model neurons (bottom; "Model").Black squares highlight broad RGC response types according to (6): OFF cells, (G1,5) ON-OFF cells (G10), fast ON cells (G18,20), slow ON (G21,23,24) and ON contrast suppressed (G28) cells, and OFF suppressed cells (G31,32).

Figure 5 .
Figure 5. Electrical single-cell recordings of responses to MEI stimuli confirm chromatic selectivity of tSbC RGCs.(a) Spiking activity (top, raster plot; middle, firing rate) of a OND RGC in response to different MEI stimuli (black bar indicates MEI stimulus duration; grey rectangle marks optimisation time window, from 1 to 1.66 s).Bottom: Activation relative to mean as a function of MEI stimulus, averaged across cells (solid line, from electrical recordings, N=4; dashed line, from Ca 2+ imaging, N=11 cells).Colours as in Figure 4. (b) Like (a) but for a sustained ON α cell (G24; N=4 cells, both for electrical and Ca 2+ recordings).(c) Different ON delayed (OND/tSbC, G28) RGC (green) dye-loaded by patch pipette after cell-attached electrophysiology recording (z-projection; x-y plane).(d) Cell from (c, green) as side-projection (x-z), showing dendritic stratification pattern relative to choline-acetyltransferase (ChAT) amacrine cells (tdTomato, red) within the inner plexiform layer (IPL).

Figure 6 .
Figure 6.Chromatic contrast selectivity of G28 RGCs derives from a nonlinear transformation of stimulus space (a) Distribution of green and UV MEI centre contrast for a linear-nonlinear (LN) model (red) and a nonlinear CNN model (black).Colour-opponent cells highlighted by filled marker.(b,c) Left: MEIs for an example cell of RGC group G28, generated with the LN model (b) or the CNN model (c).The cell's MEI centre contrast for both models is marked in (a) by asterisks.Right: Respective tuning maps of example model neuron in chromatic contrast space.Colours represent responses in % of maximum response; arrows indicate the direction of the response gradient across chromatic contrast space.(d) Difference in response predicted between LN and CNN model (in % of maximum response).(e) Contour plot of activity vs. green and UV contrast for an example tSbC (G28) RGC measured in whole-cell current-clamp mode.Labels on the contour plot indicate spike count along isoresponse curves.(f) Traces are examples of responses at the 8 extremes of -100%, 0, or 100% contrast in each colour channel.

Figure 7 .
Figure 7. Chromatic contrast tuning allows detection of ground-to-sky transitions (a) Distribution of green and UV contrasts of all movie inter-clip transitions (centre), separately for the 4 transition types, for each of which an example is shown: ground-to-sky (N=525, top left, red triangle), ground-to-ground (N=494, top right, green disk), sky-to-ground (N=480, bottom left, black downward triangle), and sky-to-sky (N=499, bottom right, purple square).Images show last and first frame of pre-and post-transition clip, respectively.Traces show mean full-field luminance of green and UV channels in last and first 1 s of pre-and post-transition clip.Black trace shows luminance averaged across colour channels.(b) Distributions as in (a), but shown as contours indicating isodensity lines of inter-clip transitions in chromatic contrast space.Density of inter-clip transitions was estimated separately for each type of transition from histograms within 10 × 10 bins that were equally spaced within the coloured boxes.Four levels of isodensity for each transition type shown, with density levels at 20 % (outermost contour, strongest saturation), 40 %, 60 % and 80 % (innermost contour, weakest saturation) of the maximum density observed per transition: 28 sky-to-ground (black), 75 ground-to-ground (green), 42 sky-to-sky (purple) and 45 ground-to-sky (red) transitions per bin.Orange markers indicate locations of N=36 G28 MEIs in chromatic contrast space (cf.Figure 3i).(c) Tuning map of G28 RGCs (N=78), created by averaging the tuning maps of the individual RGCs, overlaid with outermost contour lines from (b) (cf. Figure 6-figure supplement 2VIIb).(d,e) Same as (c) for G21 ((g), N=97) and G 5 ((h), N=33).(f) Top: Illustration of ROC analysis for two RGCs, a G21 (left) and a G28 (right).For each RGC, responses to all inter-clip transitions were binned, separately for ground-to-sky (red) and all other transitions (grey).Middle: Sliding a threshold d across the response range, classifying all transitions with response > d as ground-to-sky, and registering the false-positive-rate (FPR) and true-positive-rate (TPR) for each threshold yields an ROC curve.Numbers in brackets indicate (FPR, TPR) at the threshold indicated by vertical line in histograms.Bottom: Performance for each cell, quantified as area under the ROC curve (AUC), plotted as distribution across AUC values for all cells (black), G21 (grey), G5 (blue), and G28 (orange); AUC mean ± SD indicated as dots and horizontal lines above histograms.(g) Boxplot of AUC distributions per cell type.The box extends from the first quartile (Q1) to the third quartile (Q3) of the data; the line within a box indicates the median.The whiskers extend to the most extreme points still within [Q1 − 1.5 × IQR, Q3 + 1.5 × IQR], IQR = inter-quartile range.Diamonds indicate points outside this range.All elements of the plot (upper and lower boundaries of the box, median line, whiskers, diamonds) correspond to actual observations in the data.Numbers of RGCs for each type are indicated in the plot.(h) Illustration of stimulus with transitions as in (a) but at different velocities (50, 150, 250, and 350°/s).(i) Like (g) but for model cells and transition movies from (h) at 50°/s.(j) AUC as function of transition velocity for example RGC groups (G (1,5), (10), (18,20), (21, 23, 24), (28, 31, 32) ).
figure supplement 6VIIb), and then averaged the resulting single-cell tuning maps for each RGC group (for examples, see Figure7c-e).G 28 is most strongly tuned to fullfield transitions in the upper left quadrant containing mostly ground-to-sky inter-clip transitions (Figure7c) -unlike, for example, non-colour-opponent reference RGC groups from the slow ON and OFF response regime (Figure7d,e).Could a downstream visual area detect ground-to-sky visual context changes based on input from G 28 RGCs?To answer this question, we performed a linear detection analysis for each RGC by sliding a threshold across its responses to the inter-clip transitions, classifying all transitions that elicited an above-threshold response as ground-to-sky, and evaluating false-positive and true-positive rates (FPR and TPR, respectively) for each threshold (Figure7f).Plotting the resulting TPRs for all thresholds as a function of FPRs yields a receiver operating characteristic (ROC) curve(47) (Figure7f, middle).The area under this curve (AUC) can

800(
ACSF) solution containing (in mM): 125 NaCl, 2.5 KCl, cose, and 0.5 L-glutamine at pH 7.4.Next, the retinae were bulk-electroporated with the fluorescent Ca 2+ indicator Oregon-Green BAPTA-1 (OGB-1), as described earlier(73).In brief, the dissected retina was flat-mounted onto an Anodisc (#13, 0.2 µm pore size, GE Healthcare) with the RGCs facing up, and placed between a pair of 4-mm horizontal plate electrodes (CUY700P4E/L, Nepagene/Xceltis).A 10µl drop of 5 mM OGB-1 (hexapotassium salt; Life Technologies) in ACSF was suspended from the upper electrode and lowered onto the retina.Next, nine pulses (≈ 9.2 V, 100 ms pulse width, at 1 Hz) from a pulse generator/wideband amplifier combination (TGP110 and WA301, Thurlby handar/Farnell) were applied.Finally, the tissue was placed into the microscope's recording chamber, where it was perfused with carboxygenated ACSF (at ≈ 36°C) and left to recover for ≥ 30 min before recordings started.To visualise vessels and damaged cells in the red fluorescence channel, the ACSF contained ≈ 0.1 µM Sulforhodamine-101 (SR101, Invitrogen)(74).All procedures were carried out under dim red (> 650 nm) light.For electrophysiology experiments, we used ChAT-Cre (JAX 006410) x Ai14 (JAX 007914) mice on a C57Bl/6J background (n=2, male, aged 27 and 30 weeks).Mice were housed with siblings in groups up to 4, fed normal mouse chow and maintained on a 12:12 h light/dark cycle.Before the experiment, mice were dark-adapted overnight and sacrificed by cervical dislocation.Retinal tissue was isolated under infrared illumination (900 nm) with the aid of night-vision goggles and IR dissection scope attachments (BE Meyers).Retinal orientation was identified using scleral landmarks(75), and preserved using relieving cuts in cardinal directions, with the largest cut at the dorsal retina.Retinas were mounted on 12mm poly-D-lysine coated glass affixed to a recording dish with grease, with the GCL up.Oxygenation was maintained by superfusing the dish with carboxygenated Ames medium (US Biological, A1372-25) warmed to 32 °C.For cell-attached single cell recordings, we used Symphony software (https: //symphony-das.github.io/)with custom extensions (https:// github.com/Schwartz-AlaLaurila-Labs/sa-labs-extension).Owing to the exploratory nature of our study, we did not use randomisation and blinding.No statistical methods were used to predetermine sample size.Two-photon calcium imaging.We used a MOM-type two-photon microscope (designed by W. Denk; purchased from Sutter Instruments) (74, 76), which was equipped with a mode-locked Ti:Sapphire laser (MaiTai-HP DeepSee, Newport Spectra-Physics) tuned to 927 nm, two fluorescence detection channels for OGB-1 (HQ 510/84, AHF/Chroma) and SR101 (HQ 630/60, AHF), and a water immersion objective (CF175 LW D × 16/0.8W,DIC N2, Nikon, Germany).Image acquisition was performed with custom-made software (ScanM by M. Müller and T.E.) running under IGOR Pro 6.3 for Windows (Wavemetrics), taking time-lapsed 64 × 64 pixel image scans (≈ (100 µm) 2 at 7.8125 Hz (Figure 1c).For simplicity, we refer to such a time-lapsed scan of a local population of GCL cells as a "recording".Despite the low frame rate, the Ca 2+ re-860 sponses can be related to the spike rate (77-80).For doc-861 umenting the position of the recording fields, the retina un-862 der the microscope was oriented such that the most ventral 863 edge pointed always towards the experimenter.In addition, 864 higher resolution images (512 × 512 pixel) were acquired 865 and recording field positions relative to the optic nerve were 866 routinely logged.867 Data preprocessing.Ca 2+ traces were extracted for in-868 dividual ROIs as described previously (6, 20).Extracted 869 traces c raw were then detrended to remove slow drifts in the 870 recorded signal that were unrelated to changes in the neural 871 response.First, a smoothed version of the traces, c smooth , 872 was calculated by applying a Savitzky-Golay filter of 3 rd 873 polynomial order and a window length of 60 s using the 874 SciPy implementation scipy.signal.savgol_fil-875ter.This smoothed version was then subtracted from the 876 raw traces to yield the detrended traces.877c detrend = c raw − c smoothTo make traces non-negative (c + ), we then clipped all 878 values smaller than the 2.5 th percentile, η 2.5 , to that value, 879 and then subtracted η 2.5 from the detrended traces:880 c + = c detrend − η 2.5This procedure (i.e.clipping to, and subtracting η 2.5 ) was 881 more robust than simply subtracting the minimum.882 Finally, traces were then divided by the standard devi-883 ation within the time window before stimulus start at t 0 : 884 c := c f inal = c nn SD(c + [:t 0 ] ) For training the model on movie response, we then es-885 timated firing rates r from the detrended Ca 2+ traces c using 886 the package C2S (https://github.com/lucastheis/c2s,Theis 887 et al. (80)).888 Inclusion criteria.We applied a sequence of quality filter-889 ing steps to recorded cells before analysis illustrated in Fig-890 ure 1-figure supplement 1Ic.As a first step, we applied a 891 general response quality criterion, defined as a sufficiently 892 reliable response to the Moving bar stimulus (as quantified 893 by a quality index QI M B > 0.6), or a sufficiently reliable 894 response to the chirp stimulus (as quantified by a quality 895 index QI chirp > 0.35).The quality index is defined as in 896 ref.(6): 897 QI = Var[ r i ] t Var[r] t i where r is the T by I response matrix (time samples 898 by stimulus repetitions) and x and Var[] x denote the mean 899 and variance across the indicated dimension x, respectively.900 The second and third step made sure only cells were 901 included that were assigned to a ganglion cell group (i.e., 902 group index between 1 and 32) with sufficient confidence.903 Höfling et al. | A chromatic feature detector in the retina signals visual context changes bioRχiv | 13 to their green and UV spatial component (see Methods section Concentric anisotropic 2D Difference-of-Gaussians fit).For the analysis of MEI properties (temporal frequency, centre size, chromatic contrast), we only included cells with a sufficient DoG goodness-of-fit, determined as a value of the cost function of < .11for both green and UV on the resulting DoG fit.This threshold was determined by visual inspection of the DoG fits and led to the inclusion of 1613 out of 1947 RGCs in the MEI property analysis.Visual stimulation.For light stimulation (imaging experiments), we projected the image generated by a digital light processing (DLP) projector (lightcrafter DPM-FE4500MKIIF, EKB Technologies Ltd) through the objective onto the tissue.The lightcrafter featured a light-guide port to couple in external, band-pass filtered UV and green LEDs (light-emitting diodes) (green: 576 BP 10, F37-576; UV: 387 BP 11, F39-387; both AHF/Chroma) (82).To optimise spectral separation of mouse M-and S-opsins, LEDs were band-pass filtered (390/576 dual-band, F59-003, AHF/Chroma).LEDs were synchronised with the microscope's scan retrace.Stimulator intensity (as photoisomerization rate, 10 3 P * s −1 per cone) was calibrated to range from ≈ 0.5 (black image) to ≈ 20 for M-and S-opsins, respectively.Additionally, we estimated a steady illumination component of ≈ 10 4 P * s −1 per cone to be present during the recordings because of two photon excitation of photopigments (74, 76).Before data acquisition, the retina was adapted to the light stimulation by presenting a binary noise stimulus (20 × 15 matrix, (40 µm) 2 pixels, balanced random sequence) at 5 Hz for 5 min to the tissue.

990A 72 ×16 × 5 × 5 ×
mouse cam movie frame contained a circular field 991 of view (FOV) of 180°corresponding to 437 pixels along 992 the diameter.To minimise the influence of potential chro-993 matic and spatial aberrations introduced by the lenses, we 994 focused on image cut-outs (crops; 30°× 26°, equivalent to 995 64 pixels in size) from upper and lower visual field, 996 centred at [28°, 56°] and [−42°, −31°], respectively, rela-997 tive to the horizon (for details, see (5)).Our stimulus movie 998 consisted of 113 movie clips, each 150 frames (= 5 s) long.999 108 clips were randomly reordered for each recording and 1000 split into two 54 clips-long training sequences.The remain-1001 ing 5 clips formed a fixed test sequence that was presented 1002 before, in between, and after the training sequences (Fig-1003 ure 1b).To keep intensity changes at clip transitions small, 1004 we only used clips with mean intensities between 0.04 and 1005 0.22 (for intensities in [0, 1]).For display during the experi-1006 ments, intensities were then mapped to the range covered by 1007 the stimulator, i.e. [0, 255].1008 Convolutional neural network model of the retina.We 1009 trained a convolutional neural network (CNN) model to pre-1010 dict responses of RGCs to a dichromatic natural movie.The 1011 CNN model consisted of two modules, a convolutional core 1012 that was shared between all neurons, and a readout that was 1013 specific for each neuron (85).1014 tional neural network with 16 feature channels in each layer.1016 Both layers consisted of space-time separable 3D convolu-1017 tional kernels followed by a batch normalisation layer and 1018 an ELU (exponential linear unit) nonlinearity.In the first 1019 layer, sixteen 2 × 11 × 11 × 21 (c=#input channels (green 1020 and UV) × h=height × w=width × t=#frames) kernels were 1021 applied as valid convolution; in the second layer, sixteen 1022 11 kernels were applied with zero padding 1023 along the spatial dimensions.We parameterised the tempo-1024 ral kernels as Fourier series and added one time stretching 1025 parameter per recording to account for inter-experimental 1026 variability affecting the speed of retinal processing.More 1027 precisely, every temporal kernel was represented by the first 1028 k sine and cosine functions, with trainable weights and 1029 phases, on an evenly spaced temporal grid, where k = 7 for 1030 the first layer, and k = 3 for the second layer.Addition-1031 ally, we introduced a trainable stretch parameter for every 1032 recording to account for faster and slower response kernels.1033 For example, the first layer temporal kernels are 21 steps 1034 long.Then, in order to stay well under the Nyquist limit, 1035 we parameterise the kernels with k = 21/3 = 7 sines and 1036 cosines.1037 For each of those sines and cosines a weight (α, β) is 1038 learned to represent the shape of the temporal responses ker-1039 nel (shared among cells within a recording).Per scan i, the 1040 time grid t (21 steps from 0 to 1) is stretched by a factor τ i to 1041 account for different response speeds.To avoid adding ad-1042 ditional cycles (e.g., for stretch factors τ > 1) this is masked −(t + 21•0.95τ )

1048
In the readout, we modelled each cell's spatial recep-1049 tive field (RF) as a 2D isotropic Gaussian, parameterised as 1050 N (µ x , µ y ; σ).We then modelled the neural response as an 1051 affine function of the core feature maps weighted by the spa-1052 tial RF, followed by a softplus nonlinearity.1053 For the linearised version of the model, the architec-1054 ture was exactly the same except for the fact that there was 1055 no ELU nonlinearity after both convolutional layers.The 1056 resulting CNN was therefore equivalent to an LN model.1057 Model training and evaluation.We trained our network 1058 by minimising the Poisson loss 1059 N n=1 r(n) − r (n) log r(n)

1060
sured and r(n) the predicted firing rate of neuron n for an 1061 input of duration t=50 frames.We followed the training 1062 schedule of Lurz et al. (87).Specifically, we used early stop-1063 ping (88) on the correlation between predicted and measured 1064 neuronal responses on the validation set, which consisted of 1065 15 out of the 108 movie clips.If the correlation failed to 1066 increase during any 5 consecutive passes through the entire 1067 training set (epochs), we stopped the training and restored 1068 the model to the best performing model over the course of 1069 training.We went through 4 cycles of early stopping, restor-1070 ing the model to the best performing, and continuing train-1071 ing, each time reducing the initial learning rate of 0.01 by a 1072 learning rate decay factor of 0.3.Network parameters were 1073 iteratively optimised via stochastic gradient descent using 1074 the Adam optimiser (89) with a batch size of 32 and a chunk 1075 size (number of frames for each element in the batch) of 50.1076 For all analyses and MEI generation, we used an ensemble 1077 of models as described in ref. (34).Briefly, we trained 5 in-1078 stances of the same model initialised with different random 1079 seeds.Inputs to the ensemble model were passed to each 1080 member and the final ensemble model prediction was ob-1081 tained by averaging the outputs of the 5 members.For ease 1082 of notation, we thus redefine r(n) to be the ensemble model 1083 prediction.1084 After training, we evaluated model performance for 1085 each modelled neuron n as the correlation to the mean, i.e. 1086 the correlation between predicted response r(n) and mea-1087 sured response r (n) to the held-out test sequence, the latter 1088 averaged across 3 repetitions i = {1, 2, 3}: C(r (n) , r (n) i i ).1089 Unlike the single-trial correlation C(r (n) , r (n) i ) which is al-1090 ways limited to values < 1 by inherent neuronal noise, a per-1091 fect model can in theory achieve a value of 1 for the corre-1092 lation to the mean, in the limit of infinitely many repetitions 1093 when the sample average r (n) i i is a perfect estimate of the 1094 true underlying response ρ (n) .The observed correlation to 1095 the mean can thus be interpreted as an estimate of the frac-1096 tion of the maximally achievable correlation achieved by our 1097 model.For deciding which cells to exclude from analysis, 1098 we used average single-trial correlation ( C(r (n) , r (n) i ) i ) 1099 since this measure reflects both model performance as well 1100 as reliability of the neuronal response to the movie stimulus 1101 for neuron n (see also Methods section on Inclusion crite-1102 ria).1103 Synthesising MEIs.We synthesised maximally exciting 1104 inputs for RGCs as described previously (32).Formally, for 1105 each model neuron n we wanted to find 1106 x * (n) = arg max x r(n) (x) 30:50 t , (3) i.e. the input x * (n) where the model neuron's re-1107 sponse r(x) 30:50 t , averaged across frames 30 to 50, at-1108 tains a maximum, subject to norm and range constraints 1109 (see below).To this end, we randomly initialised an input x (n) 0 ∈ R c×w×h×t of duration t=50 frames with Gaussian 1111 Höfling et al. | A chromatic feature detector in the retina signals visual context changes bioRχiv | 15 colour channel separately and ensured that the range of the MEI values stayed within the range covered by the training movie.This was achieved by clipping values of the MEI exceeding the range covered by the training movie to the minimum or maximum value.Optimisation was run for at least 100 iterations, and then stopped when the number of iterations reached 1,000, or when it had converged (whichever occurred first).Convergence was defined as 10 consecutive iterations with a change in model neuron activation of less than 0.001; model neuron activations ranged from ≈ 1 to ≈ 10.We denote the resulting MEI for neuron n as x * (n) .

Spatial and temporal components
of MEIs.For each colour channel c, we decomposed the spatiotemporal MEIs into a spatial component and a temporal component by singular value decomposition: U, S, V = svd(x * (n) c ) with x * (n) c ∈ R 50×288 for c ∈ [green, UV] is the MEI of neuron n in a given colour channel with its spatial dimension (18×16 = 288 ) flattened out.As a result, any spatiotemporal dependencies are removed and we only analyse spatial and temporal properties separately.The following procedures were carried out in the same manner for the green and the UV component of the MEI, and we drop the colour channel index c for ease of notation.The temporal component is then defined as the first left singular vector, U :1 , and the spatial component is defined as the first right singular vector, V T :1 , reshaped to the original dimensions 18 × 16.Concentric anisotropic 2D Difference-of-Gaussians fit.We modelled the spatial component as concentric anisotropic Difference-of-Gaussians (DoG) using the nonlinear leastsquares solver scipy.optimize.least_squareswith soft-L1 loss function (40).The DoGs were parameter-1156 ized by a location (µ x , µ y ) shared between centre and sur-1157 round, amplitudes A c , A s , variances (σ c x , σ c y ), (σ s x , σ s y ), and 1158 rotation angles θ c , θ s separately for centre and surround:

1162
lowing way: Since we set the model readout's location pa-1163 rameters to (0, 0) for all model neurons when generating 1164 their MEIs, we also expected the MEIs to be centred at (0, 1165 0), as well.Hence, we determined the location of the min-1166 imum and the maximum value of the MEI; whichever was 1167 closer to the centre (0,0) provided the initial values for the 1168 parameters (µ x , µ y ).Starting from there, we then first fit a 1169 single Gaussian to the MEI, and took the resulting param-1170 eters as initial parameters for the DoG fit.This was a con-1171 strained optimisation problem, with lower and upper bounds 1172 on all parameters; in particular, such that the location param-1173 eter would not exceed the canvas of the MEI, and such that 1174 the variance would be strictly positive.1175 MEI properties.1176 Centre size We defined the diameter of the centre of 1177 the MEI in the horizontal and the vertical orientation, 1178 respectively, as d c x = 2σ c x and d c y = 2σ c y .The centre size 1179 was calculated as 1 2 (d c x + d c y ).We then estimated a contour 1180 outlining the MEI centre as the line that is defined by all 1181 points at which the 2D centre Gaussian G c attains the value 1182 G c (x, y) with (x, y) = (µ x + σ c x , µ y + σ c y ).The centre mask 1183 m was then defined as a binary matrix with all pixels within 1184 the convex hull of this contour being 1 and all other pixels 1185 set to 0. This mask is used for calculating centre chromatic 1186 contrast (see below).

1187 1188
Temporal frequency To estimate temporal frequency of 1189 the MEIs, we estimated the power spectrum of the temporal 1190 components using a Fast Fourier Transform after attenuating 1191 high frequency noise by filtering with a 5 th order low-pass 1192

( 1 ,
were randomly interleaved.With 11 stimuli, presented at 25 positions and lasting 2 s each, the total stimulus duration was 11 × 25 × 2 s = 550 s.Since the model operated on a z-scored (0 mean, 1 SD) version of the movie, MEIs as predicted by the model lived in the same space and had to be transformed back to the stimulator range ([0, 255]) before being used as stimuli in an experiment by scaling with the movie's SD and adding the movie's mean.The MEIs' green channel was then displayed with the green LED, and the UV channel was displayed with the UV LED.For experiments at Northwestern University, an additional transform was necessary to achieve the same levels of photoreceptor activation (photoisomerization rates) for M-and S-cones with different LEDs.To ensure proper chromatic scaling between the different experimental apparatuses with different spectral profiles, we described the relative activation of M-and S-cones by the green and UV LEDs in the stimulation setup used in the two photon imaging experiments (setup A) by a matrix A = a mg a sg a mu a su = and the relative activation of M-and S-cones by 1241 the stimulation setup used in the patch-clamp experiments 1242 where diagonal entries describe the activation of M-1244 cones by the green LED, and of S-cones by the UV LED, 1245 and entries in the off-diagonal describe the cross-activation 1246 (i.e., M-cones by UV-LED and S-cones by green LED).The 1247 activation of M-cones and S-cones e T = (e m , e s ) by a stim-1248 ulus x ∈ R 2×1 displayed on a given stimulation setup was 1249 approximated as e = Ax (90).Hence, a stimulus x dis-1250 played on setup B, defined as x = B −1 Ax, will achieve 1251 the same photoreceptor activation as stimulus x displayed 1252 on setup A. Since the solution exceeded the valid range of 1253 the stimulator ([0, 255]), we added an offset and multiplied 1254 with a scalar factor to ensure all stimuli were within the valid 1255 range.1256 Analysing RGC responses to MEI stimuli.We wanted to 1257 evaluate the responses of RGCs to the MEI stimuli in a spa-1258 tially resolved fashion, i.e. weighting responses to MEIs 1259 displayed at different locations proportional to the strength 1260 of the RGCs RF at that location.In order to be able to 1261 meaningfully compare MEI responses between RGCs and 1262 across groups, for each RGC, we first centred and scaled 1263 the responses to zero mean and a standard deviation of 1. 1264 Then, for each RGC n, we computed a spatial average of its 1265 responses, weighting its responses at each spatial location 1266 (x, y) proportional to the Gaussian density N µ n ,σ n (x, y), 1267 where the parameters of the Gaussian µ n = (µ x , µ y ), σ n 1268 were the model's estimated readout parameters for neuron 1269 n (Figure 4b,c,d left): 11×60 is the 60 frames (2 s) long re-1271 sponse of neuron n to an MEI at position (x, y) = (x , y ), 1272 resampled from the recording frame rate of 7.81 Hz to 1273 30 Hz.We then averaged r (i) x,y across time in the op-1274 timisation time window, i.e. frames 30-50, to get a scalar 1275 response r(n) = r (n) x,y,t for each MEI stimulus (Fig-1276 ure 4d).

(
1284 difference (in units of SD response) between the response to 1285 the MEI of interest (x * i ) and the mean response to all other 1286 r(n) (x * j ), averaged across all cells 1287 n belonging to the group of interest G g .1288 Characterising nonlinear processing of chromatic 1289 contrast space.We wanted to analyse the tuning of 1290 G 28 /tSbC RGCs to chromatic contrast and to this end, we 1291 mapped the model response and its gradient across chro-1292 matic contrast space (Figure 6).Specifically, the MEIs 1293 have d = 2 × 18 × 16 × 50 = 28, 800 pixels and dimensions, 1294 14, 400 for each colour channel.Now let x * (n) ∈ R 1x28800 1295 be the cell's MEI estimated using the LN model, with the 1296 first d=14,400 dimensions defining the green pixels and the 1297 remaining dimensions defining the UV pixels.Then for each 1298 cell we consider a two-dimensional subspace spanned by 1299 two basis vectors e 1 , e 2 where 1300

1301
figure supplement 3VIIIa,b).Each edge of the trajectory is 1353

1395withFigure I 1 -
Figure I 1-figure supplement 1.(a) Distribution across cell types for this dataset, and for the dataset described in Baden et al. (6) which was the basis for our classifier (81).(b) Mean ± SD of model performance, evaluated as correlation between model prediction and RGC response on the 25 s long test sequence, averaged across 3 repetitions of the test sequence, for each cell type.(c) Response quality, RGC group assignment and model performance filtering pipeline showing the consecutive steps and the fraction of cells remaining after each step.Black bars and numbers indicate cells from all experiments (i.e.all RGCs for which we recorded chirp, MB, and movie responses), red bars and numbers indicate the subset of cells recorded in the MEI validation experiments (i.e.those RGCs for which we additionally recorded MEI stimuli responses).Dotted bars indicate the number of cells before the current filtering step.The filtering steps were as follows: 1. Keep only cells that pass the chirp OR MB quality criterion (QI M B > .6OR QI chirp > .35). 2. Keep only cells that the classifier assigns to a group with confidence ≥ .25. 3. Keep only cells assigned to a ganglion cell group (groups 1-32; groups 33-46 are amacrine cell groups); 4. Keep only cells with sufficiently high model performance ( C(r (n) , r (n) i

Figure II 3-figure supplement 1 .
Figure II 3-figure supplement 1. Example MEIs for example cell types.Rows in each panel as in Figure 4a.

Figure IV 4 - 25 Figure V 5 -
Figure IV 4-figure supplement 1.(a) Recorded (top, r) and predicted (bottom, r) responses to the 11 different MEIs for all example cell types.Left: responses are averaged across the indicated dimensions x, y, n: different MEI locations (x, y) and RGCs in a group (n); black bar indicates stimulus duration (from 0 to 1.66 s), grey rectangle marks optimisation time window (from 1 to 1.66 s).Right: Responses to different MEIs, additionally averaged across time (t) within the optimisation time window.Error bars indicated SD across cells.(b) Correlation between the measured and predicted response magnitudes to the MEI stimuli per example cell type.Cumulative histogram is across all N=788 cells; 50% of cells have a correlation between measured and predicted response magnitude of ≥ 0.8.(c) Mean ± SD of selectivity index (see Methods) for the example cell groups, indicating the difference in response to MEI 28 vs. the average response to all other MEIs in units of standard deviation of the response.
(h) Boxplot of AUC distributions per cell type (dorsal).The box extends from the first quartile (Q1) to the third quartile (Q3) of the data; the line within a box indicates the median.The whiskers extend to the most extreme points still within [Q1 − 1.5 × IQR, Q3 + 1.5 × IQR], IQR = inter-quartile range.Diamonds indicate points outside this range.All elements of the plot (upper and lower boundaries of the box, median line, whiskers, diamonds) correspond to actual observations in the data.Numbers of RGCs for each type are indicated in the plot.

Figure VIII 6 -
Figure VIII 6-figure supplement 3. (a) Illustration transition stimulus paradigm (from Figure 7h).(b) Structure of stimuli for different velocities, using a ground-to-sky transition as an example.(c) Statistics of the area under the ROC curve (AUC) for the sky-ground detection task in the simulation for different velocities (G28 vs. the next-best RGC group).Columns (from left): mean ± standard deviation of AUC values (top: G28; bottom: the respective best next RGC type); difference in mean AUC and corresponding bootstrapped 95% confidence intervals; Cohen's d and p-value of a two-sample permutation test with 100,000 repeats.(d) Boxplots of AUC distributions per cell type for the different velocities (plots like in Figure 7g,j).

representation allows for detection of 461 ground-to-sky transitions. Studies analysing visual
460 Warped 462 scenery from the mouse's perspective have repeatedly 463 found that chromatic contrast changes strongly at the 464 We wanted to test how likely the differ-, φ ν g ) for each sample.We then estimated the 95 1392 % confidence interval for Γ(φ s , φ ν g ) as the interval defined 1393 by the 2.5 th and 97.5 th percentile of the bootstrapped dis- ). 1365 Statistical analysis.1366 Permutation test.1372 ∆AUC, by shuffling the RGC group labels of the two groups 1373 of interest (e.g.G 28 and G 24 ) and calculating the test statis-1374 tic with shuffled labels 100,000 times.We only included 1375 RGC groups with at least N=4 cells in this analysis.We then 1376 obtained a p-value for ∆AUC observed with true labels as 1377 the proportion of entries in the null distribution larger than 1378 ∆AUC.1379 Bootstrapped confidence intervals.We bootstrapped con-1394 tribution of Γ(φ s , φ ν g ).